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STRINGS IN THE 77 CARINAE NEBULA: HYPERSONIC RADIATIVE COSMIC BULLETS 
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ABSTRACT 

We present the results of a numerical study focusing on the propagation of a hypersonic bullet subject to 
radiative cooling. Our goal is to explore the feasibility of such a model for the formation of "strings" observed 
in the -q Carinae Homunculus nebula. Our simulations were performed in cylindrical symmetry with the adaptive 
mesh refinement code AstroBEAR. The radiative cooling of the system was followed using the cooling curve by 
Dalgarno & McCray (1972). In this letter we discuss the evolution and overall morphology of the system as well 
as key kinematic properties. We find that radiative bullets can produce structures with properties similar to those 
of the 77 Carinae strings, i.e. high length-to-width ratios and Hubble-type flows in the form of a linear velocity 
increase from the base of the wake to the bullet head. These features, along with the appearance of periodic "ring- 
like" structures, may also make this model applicable to other astrophysical systems such as planetary nebulae, 
e.g. CRL 618 and NGC 6543, young stellar objects, etc. 

Subject headings: circumstellar matter — stars: individual (r; Carinae) — stars: mass-loss — ISM: jets and 
outflows — planetary nebulae: individual (CRL 618) 



1. INTRODUCTION 

Of the many intriguing questions surrounding the luminous 
blue variable 77 Carinae, the origin of the long thin "strings" 
found in the outer nebular regions remains particularly vexing. 
First observed by Meaburn et al. (1996) and further studied by 
Weis, Duschl, & Chu (1999) the strings exhibit two remarkable 
properties: very high values of the length-to-width ratios and 
the presence of the Hubble-type flow, i.e. linear velocity in- 
crease from the base of the strings to their tip. The strings also 
show a surface brightness decrease toward the string head and 
the true tips possibly may be invisible in optical wavelengths. 
These features represent the key challenges that must be met by 
any model that hopes to offer a successful explanation of the 
strings' origin and evolution. 

Several models were suggested for the strings including jets 
(Garcfa-Segura et al. 1999) and ionization shadows (Soker 
2001), see Redman, Meaburn, & Holloway (2002) for a com- 
plete discussion of the suggested models. A particularly 
promising model has been proposed by Redman et al. (2002) 
who suggest that a single hypersonic bullet propagating in the 
ambient medium is generating the strings observed in the r\ 
Carinae nebula. Since the Great Eruption that ejected the main 
nebula was know to be an impulsive event with most of the mass 
and momentum directed in the polar directions (Smith et al. 
2003), and since the strings seem to be coeval with the overall 
nebula, it is plausible that fragments or bullets were produced 
in that same event. 

In this work we do not address the issue of bullet origin fo- 
cusing instead on the nature of their evolution and observable 
properties. In particular, we investigate the feasibility of the 
Redman et al. (2002) model and we address two questions: (1) 
can a radiatively cooled hypersonic bullet produce structures 
with large length-to-width ratios, and (2) do such systems ex- 
hibit the Hubble-type flow behaviour in their downstream wake. 

In Section 2 we briefly discuss the numerical code used, the 
key dimensionless parameters characterizing the system, and 



the setup of the numerical simulation that was carried out. In 
Section 3 we discuss the results of our numerical study, and 
in Section 4 we present the conclusions including possible an- 
swers to the questions posed above as well as applicability of 
our results to other astrophysical systems. 

2. NUMERICAL STUDY 

The simulation presented in this paper was performed with 
the AstroBEAR code based on the BEARCLAW adaptive mesh 
refinement package (Mitran 2001; Berger & LeVeque 1998). 
The code allows for arbitrary levels of refinement within com- 
puter memory constraints, as well as adaptive numerics and 
multiphysics. The Euler hydrodynamic equations with cooling 
source terms were solved using explicit second-order-accurate 
scheme. The wave propagation method (LeVeque 1997) was 
applied with the full nonlinear Riemann solver. An oper- 
ator split set of ODE's, describing the source term effects, 
was solved using the fourth-order-accurate fully implicit Kaps- 
Rentrop-type scheme (Kaps & Rentrop 1979) (see also (Press 
etal. 1997)). 

There are three dimensionless parameters that describe the 
evolution of the system studied. The first two are the den- 
sity contrast between the bullet and the ambient medium \b = 
Pb I Pa and the Mach number Mb of the bullet velocity at time 
t = 0.0. Those two parameters determine the overall hydro- 
dynamic regime of the system evolution. The third one, the 
cooling parameter ipt,, describes the effects of cooling and is 
defined as 

A=^- (1) 

thydro 

Here the hydrodynamic timescale t^dro is defined as the bul- 
let crushing time (or clump crushing time t cc as described by 
Poludnenko, Frank, & Blackman (2002)) 
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where r b is bullet radius, vj is bullet velocity, F c \ « 1.3 and 

2.16 

^ 1 + IT6^P72- (3) 

The two factors F c i and F sr relate the unperturbed upstream 
conditions with the internal bullet post-shock ones and are de- 
scribed in (Poludnenko et al. 2002). The interaction of the bul- 
let with the ambient medium launches a shock propagating in- 
side the bullet with the post-shock temperature T ps and p ps . The 
cooling rate A ps for such a temperature can be found based on 
the cooling curve given by Dalgarno & McCray (1972). Thus 
the cooling timescale is estimated as follows 

tcool = T ' ,s Tmm fc OTg; (4) 

where T min = 100 K is the minimum temperature gas can cool 
down to. The internal bullet post-shock temperature T ps and 
density p ps can be found in the usual manner using the Rankine- 
Hugoniot relations. T ps and density p ps are the functions of the 
unshocked bullet temperature and density respectively as well 
as of the internal bullet shock Mach number M, s described by 
the expression 

M is =M b (F cl F st ) l/2 . (5) 

We carried out simulations covering bullet Mach numbers in 
the range Mb = 10-200. In terms of the cooling parameter 
we explored the regimes ranging from the quasi-adiabatic one 
ipb > 1 to the strongly cooling one tpb w 10~ 5 . In this paper we 
will report only on the case that allows us to best address the 
issues related to rj Carinae strings. From our parameter space 
explorations this turns out to be the case with the strongest cool- 
ing = 2.8 • 10~ 5 . We defer broader conclusions gleaned from 
the parameter space exploration to a subsequent paper. 

Here we present a simulation of bullet propagation where 
the initial bullet radius is r\, = 3.0 • 10 15 cm. The initial bul- 
let density is pb = 10 5 cm" 3 in the inner 80% of the bullet ra- 
dius and in the outer 20% of the radius density is smoothed out 
through the tanh-type function (cf. expression (1) in (Polud- 
nenko et al. 2002)). This translates into the initial bullet mass 
m b « 0.5 • 10" 5 M Q . The initial bullet Mach number is M b = 20, 
or vj w 235 km s" 1 . The computational domain is of size 
1.8 • 10 17 cm x 2.4- 10 16 cm or expressed in terms of bullet 
radii 60 r b x 8r&. The ambient density is p a = 1000 cm" 3 and 
the ambient temperature T a = \Q A K. It is assumed that initially 
the bullet is in pressure equilibrium with the ambient material. 

The run was carried out in cylindrical symmetry with the x- 
axis being the symmetry axis. We use an adaptive grid with 4 
levels of refinement. This provided a resolution of 128 cells per 
cloud radius or an equivalent resolution of 1024 x 7680 cells. 
Although the above resolution is sufficient to place the simula- 
tion in the converged regime in the adiabatic case (Klein, Mc- 
Kee, & Colella 1994), in the cooling case one has to be more 
precise in one's definition of convergence. Since, as it will be 
discussed below, the key process in the system is bullet frag- 
mentation via instabilities at the upstream surface, our criterion 
of convergence was the constancy of the initial fragmentation 
spectrum. In this sense the above resolution ensured that the 
instability wavelength, and therefore the number of fragments 
into which the bullet breaks up, does not change with increas- 
ing resolution. On the other hand, it should be mentioned that 
such resolution might not be sufficient to track the contraction 

4 Note that the simulation run time is larger than the age of the rj Carinae strings 
system evolution, the results can be appropriately scaled to r) Carinae conditions. ^ 



of each fragment due to radiative cooling to its smallest size 
defined by the equilibrium between cooling and external heat- 
ing. Finally, outflow boundary conditions were prescribed on 
all boundaries. 

For the simulation described above the cooling parameter 
was tpb = 2.8 • 10" 5 . The hydrodynamic timescale, i.e. the bullet 
crushing time, was thydm = 2.54 • 10 9 s w 80.5 yrs., the cooling 
timescale was t coo i = 70.9 • 10 3 s 19.7 hours. The simulation 
total run time was 9.49 • 10 9 s w 251.1 yrs 4 . Thus the evolution 
of the system was in a regime strongly dominated by cooling. 

3. RESULTS 

3.1. Overall Morphology and Evolution of the System 

When the interaction of an inhomogeneity (bullet or clump) 
proceeds in the adiabatic regime the dominant process defining 
the evolution of the system is the lateral clump re-expansion 
(Klein et al. 1994; Poludnenko et al. 2002). The internal shock 
compresses the clump and once such compression is completed 
the re-expansion begins into regions of lowest total pressure, 
i.e. in the lateral direction. Such re-expansion, acting in com- 
bination with instabilities at the upstream clump surface, is re- 
sponsible for the destruction of the clump. When more than one 
clump is present, lateral re-expansion also drives interclump in- 
teractions via merging and formation of larger structures that 
subsequently alter the global flow (Poludnenko et al. 2002). 

The fundamental distinction between the radiatively cooled 
inhomogeneous systems, including individual bullets, and the 
adiabatic ones is the minimal role of lateral re-expansion. In- 
stead, the dominant process is the formation of instabilities at 
the upstream surface of the bullet with a wavelength signifi- 
cantly smaller than the one observed in the adiabatic case. As 
the bullet begins to drive through the ambient medium, hydro- 
dynamic instabilities (Richtmeyer-Meshkov, Rayleigh-Taylor) 
produce the initial instability seed. The resulting density vari- 
ations quickly trigger the onset of thermal instabilities which 
thereafter become the dominant process responsible for bullet 
fragmentation. In our simulation such instabilities started de- 
veloping with a wavelength of about 8 - 10% of the initial bul- 
let radius or 2.4 • 10 14 -3.0 • 10 14 cm w 16-20 a.u. This initial 
scale is crucial for subsequent evolution since it determines not 
only the continuing process of bullet fragmentation but also the 
structure of the downstream flow. As it was mentioned above, 
we did not observe any significant changes in the fragmentation 
spectrum with changing grid resolution, therefore we believe 
that we observed the true fragmentation spectrum correspond- 
ing to the given flow conditions. Of course this conclusion may 
be altered with fully 3-D simulations. The question of predict- 
ing the properties of the initial fragmentation spectrum for the 
given characteristics of the system is important for the under- 
standing of the physics of radiative hypersonic inhomogeneous 
systems, however it is outside the scope of the current paper 
and it will be investigated in the subsequent work. 

Since the wavelength of the initial fragmentation spectrum is 
approximately constant along the bullet radius, the fragments 
produced by the instability are of different mass with the most 
massive ones being closest to the symmetry axis and the out- 
ermost ones being the lightest. As a result the fragments are 
"peeled off" from the bullet one by one starting with the out- 
ermost ones in radius. Each fragmentation event results in the 

(~ 150 yrs.). However, since it is the dimensionless numbers that determine the 
Ve discuss that in further detail in Section 3.2. 
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formation of a distinct "ring-like" feature in the bullet wake. 
The formation of rings is a consequence of the axisymmetry 
of these simulations. In 3-D it is likely that the rings would 
themselves fragment (Klein et al. 2003; Robey et al. 2002). 

One noteworthy point is that the evolution of the bullet and its 
fragments suggests the presence of steady "mass loading" of the 
downstream flow via the hydrodynamic ablation. Mass loading 
has been claimed to be an important process in all clumpy hy- 
drodynamic systems (Hartquist et al. 1986; Hartquist & Dyson 
1988). The extent and nature of such mass loading will be in- 
vestigated in the subsequent work. 

Figure 1 shows the computational domain at time t = 25 1 . 1 
yrs. In Figure la the synthetic Schlieren image (gradient of the 
density logarithm) is shown illustrating the shock and vortex 
sheet structure in the flow. Figure lb shows the synthetic ob- 
servation image of the computational domain. Since our sim- 
ulation did not track the full ionization dynamics of the flow, 
the image represents the total radiative energy losses summed 
along each ray. The 2D distribution of the state vector obtained 
in the simulation was extended using cylindrical symmetry to a 
2048 x 2048 x 7680 cells uniform grid 3D data cube. Thus the 
synthetic observation image represents the 2D projected dis- 
tribution of the logarithm of the emissivity / integrated in the 
z-direction according to the formula 

(6) 



lent region immediately behind the bullet head and in between 
the stripped bullet fragments, the gas tends to re-expand essen- 
tially at the sound speed of the ambient gas and fill the cavity. 
Such re-expansion causes gas to rebound on the axis resulting 
in a reflected shock. The effect of this shock can be seen in 
Figure la around the symmetry axis as a complex shock region 
that is narrower than the bowshock. 

Note, that in Figure lb the surface brightness of the system 
increases toward the bullet head. This differs from the surface 
brightness behaviour observed in the r\ Carinae strings. As was 
noted by Weis et al. (1999) the string surface brightness de- 
creases toward the head so that the true string tip might not even 
be detected. Recall, however, that Figure lb is not a true syn- 
thetic observation image in a certain line but rather a projected 
map of total radiative energy losses in the system. In addition 
illumination by stellar photons is not included in our calcula- 
tion. Thus in the simulation the bullet head and its vicinity, 
where the density and temperature are highest so that n 2 A(T) 
are large, appear brightest. In reality gas in the bullet head and 
near it may be too hot to be visible optically as was suggested 
by Redman et al. (2002). In our further work we plan to in- 
clude tracking of full ionization dynamics that will allow us to 
produce true synthetic observation images in a particular emis- 
sion line or a combination of lines. 



where i, j, and Ic are the cell indices in the x-, y-, and z- 
direction respectively, and the cooling rate A(T) here, as well 
as in the simulation, was determined based on the cooling curve 
described by Dalgarno & McCray (1972). 

There are several ring-like structures in Figure 1 in the bul- 
let wake. These result from the fragmentation events men- 
tioned above. The most prominent ring occurs at the distance of 
w 1 .38 • 10 17 cm from the bullet head (ring 1) with the width of 
« 2.94 • 10 16 cm w 1965 a.u. which is the widest part of the bul- 
let wake. In Figure lb several other bright rings, resulting from 
fragmentation events, are visible closer to the bullet head. They 
gradually increase in radius further downstream with the radius 
of the largest ring (ring 2), located at the distance of w 5.8 • 10 16 
cm from the bullet head, being « 1.82 • 10 16 cm w 1216 a.u. 
Note that the width of the bullet wake stays practically con- 
stant, increasing only slightly, for a distance of « 8.0 • 10 16 cm 
« 5350 a.u. from ring 2 to ring 1. Thus in our simulations the 
length-to-width ratio of the bullet wake is 6.1 -9.9 and the bul- 
let wake width is 4-6.5 times larger than the width of String 
5 while they share the same length. Note also that the width 
of the bullet head is only « 1.05 • 10 15 cm w 70 a.u. We fol- 
lowed the simulation beyond the time t = 25 1 . 1 yrs, shown in 
Figures 1 and 2, up to the time t = 300.8 yrs. and we did not 
observe any appreciable changes in the aforementioned values 
of the length-to-width ratio. 

The second process that determines the structure of the flow 
in the wake is gas re-expansion into the cavity excavated by the 
bullet. The highest temperature reached by gas in the system 
is at the tip of the central bullet fragment and is about 3.5 • 10 5 
K. Gas passing through the bullet bowshock cools down very 
rapidly, weakening the bowshock and making it more oblique, 
which in its turn prevents further heating of the ambient gas. As 
it can be seen in Figure la, the bowshock practically disappears 
half-way downstream from the bullet head. Beyond the turbu- 

5 It should be noted that we observed similar behaviour of the total velocity distribution in other systems with the cooling parameter up to -0i> 
approaching the adiabatic regime with ip > 1.0 the velocity distribution takes a more complicated form. 



3.2. Velocity Distribution and Hubble-type Flow 

Another key property of the strings is the presence of 
Hubble-type flows in the bullet wake. The left panel of Figure 2 
shows the distribution of the total velocity v tot (x) along the sym- 
metry axis of the bullet/wake as a function of distance from the 
bullet head. The linear velocity decrease from the maximum 
value of 210 km s _1 from head to base is clear aside from some 
minor fluctuations arising due to the unsteadiness of the down- 
stream flow 5 . Note that very little deceleration of the bullet 
head has occurred. Given that the initial bullet velocity is 235 
km s"'we see that after 251.1 yrs., in which the bullet mate- 
rial traveled the distance of 1.8 • 10 17 cm « 0.058 pc, the bullet 
material lost only w 10% of its original velocity. That is due 
to the fact that the central fragment retains most of the original 
bullet mass while it has a rather small cross-section due to the 
cooling-induced contraction. 

Weis et al. (1999) quote a velocity slope in the r\ Cari- 
nae strings ranging from 2590 km s _1 pc _1 (String 5) to 3420 
km s _1 pc _1 (String 2). In our simulation we find a velocity slope 
of 3600 km s _1 pc _1 which is about 5% - 40% higher. Note, that 
the length of String 5 is practically equal to the distance over 
which we allowed the bullet to propagate, however String 2 is 
about 25% shorter (Weis et al. 1999). 

One should be cautious, however, in directly comparing the 
total velocity distribution shown in the left panel of Figure 2 
to the observationally determined velocity distributions in the 
strings of r\ Carinae (cf. Figure 6 in (Weis et al. 1999)). Quan- 
tities, more closely resembling the ones that are obtained obser- 
vationally, should be employed in order to make the compari- 
son between the numerical results and observations relevant. 
For example, one can create emissivity-weighed total velocity 
Vemis(x) maps of the flow. Distribution of such a quantity along 
the x-axis is shown in the right panel of Figure 2. v emis (x) was 

0.25. For systems 
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obtained according to the formula 

_ y.j>n X i v J Jy""" n 2 (x,y)A(x,y)v wt (x,y)dy 

V,M) = T^A7 M "'"J;™; nH Xl y)A( Xl y)dy ' (?) 
where / and / are the cell indices in the x- and y-directions, Vj 
is the total velocity in the cell j, and in the summation in the 
numerator we included only the cells with non-zero total veloc- 
ity. The distribution of this quantity is significantly more noisy 
than the total velocity cross-cut. However, the local maxima 
of v em i S (x) roughly tend to fall on the same line as that in the 
previous figure giving some indication of the presence of the 
Hubble-type flow. 

The distance traveled by the bullet in our simulation is prac- 
tically equal to the length of String 5 and is factors of 1 .5 - 3.0 
shorter than the rest of the strings (except for String 2) (Weis et 
al. 1999). If we assume that the true length of strings is not their 
observed length but rather the distance of the tip from the star, 
than the r\ Carinae strings are the factors of 3.2-5.56 longer 
(Weis et al. 1999) than the computational domain in our simula- 
tion. However, the kinematic age of the bullet in our simulation 
is about 67% larger than that of the r\ Carinae strings: 251.1 
years vs. 150 years. This indicates that in our simulation the 
bullet velocity was approximately a factor 2.5-5.0 (or a factor 
5.35-9.3 if we assume that the string base is located near the 
central star) smaller than in the case of r\ Carinae strings, i.e. 
the bullet velocity should be 590- 1 175 km s^or 1260-2185 
km s" 1 ) with respect to the central star. An increased bullet 
velocity would decrease the width of the bullet wake and in- 
crease the length-to-width ratio making a better match between 
our simulations and the strings of r/ Carinae. However, a sig- 
nificant increase in velocity, and therefore Mach number, re- 
sults in much higher temperatures in the bullet wake and the 
bow shock. That, in turn, increases the cooling parameter ipb 
eventually sending the system into the adiabatic regime which 
no longer exhibits the high length-to-width ratio of the wake 
and the Hubble-type downstream flow 6 . The easiest way to 
reconcile the need for the significantly higher bullet velocity, 
required in order to obtain the correct kinematic age, with the 
need for the relatively low post shock temperatures, required 
to decrease the bullet wake width, is to embed the bullet in the 
ambient flow that itself moves with significant speed relative to 
the central star. There is some observational evidence for the 
presence of such high velocity ambient flow (Weis, Duschl, & 
Bomans 2001). For example, string 1 may be embedded in the 
material moving with velocity 500-600 km s"'with respect 
to the central star, whereas the maximum velocity of the string 
material is w 995 km s -1 (\Veis et a l. 1999). 

4. CONCLUSIONS 

In this paper we have presented a numerical study of a hyper- 
sonic radiative cosmic bullet in an attempt to model the strings 



observed in the r\ Carinae Homunculus nebula. Our simulations 
follow from the discussion by Redman et al. (2002). 

Our principal conclusions are as follows: (1) hypersonic ra- 
diative bullets are capable of producing structures with high 
length-to-width ratios (between 6 and 10 for our study); (2) the 
dominant process responsible for bullet destruction is instabil- 
ity formation at the bullet upstream interface leading to sep- 
arate fragmentation episodes; these result in the formation of 
periodic "ring-like" structures in the bullet wake; (3) the sim- 
ulations do show the presence of Hubble-type flows along the 
axis in the bullet wake in terms of a linear decrease of the to- 
tal velocity downstream from the bullet head. How these flows 
appear observationally remains an open question. 

Thus we conclude that the bullet model appears as a good 
candidate for the strings of rj Carinae, though further work is 
needed. If this model finds continued success then theorists 
will confront the issue of how such high velocity bullets are 
generated, i.e. at the star or at larger distances. Confirming the 
existence of such bullets may be helpful in determining the na- 
ture of the processes, occurring within the star, which were part 
of its various eruptions. 

The astrophysical applicability of the model, discussed in 
this paper, might extend beyond the r\ Carinae strings. One 
of the most spectacular examples of nebular systems with long 
thin structures is the protoplanetary nebula CRL 618 and, in 
particular, its shocked lobes. The most prominent features of 
the lobes are the periodic rings similar to the ones present in 
our simulation. As it was discussed, such rings arise naturally 
as a consequence of the bullet fragmentation events. There is 
also some evidence for the velocity increase in the lobes from 
the base to the tip (Sanchez Contreras, Sahai, & Gil de Paz 
2002). Other examples include, but are not limited to, other 
planetary nebulae, e.g. strings in NGC 6543 (Weis et al. 1999), 
HH objects, etc. 

Further developments of the model should include the study 
of the details of the formation of the initial fragmentation spec- 
trum, investigation of the importance of mass-loading due to 
the ablation of the bullet head and bullet fragments, and inclu- 
sion of more realistic description of ionization dynamics and 
radiative cooling in the system. The latter will allow us to pro- 
duce more realistic synthetic observation images and synthetic 
spectra for better comparison with observational data. 

This work was supported in part by the NSF grant AST- 
9702484, NASA grant NAG5-8428 and the Laboratory for 
Laser Energetics under DOE sponsorship. 

The most recent results and animations of the numer- 
ical experiment, described above as well as the ones 
not mentioned in the current paper, can be found at 
www.pas.rochester.edu/~wma. 
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FIG. 1. — a) Synthetic Schlieren image of the computational domain at time 251.1 yrs. Shown is the gradient of the density logarithm, b) Synthetic observation 
image of the computational domain for the same time as in a) (see text for the detailed description). Note the periodic ring-like structures in the domain resulting 
from individual fragmentation episodes. 
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FIG. 2. — Left: distribution of total velocity along the symmetry axis of the bullet at the time 251.1 yrs. (same time in the simulation as the one shown in Figure 1). 
Right: distribution of the emissivity-weighed total velocity in the system at the same time (see text for the detailed description). 



